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1.  Introduction 


This  report  documents  a  set  of  functions,  written  in  C++,  that  can  be  used  to  generate 
pseudorandom  numbers  that  have  either  uniform  or  normal  distributions  and  pseudorandom 
integers  that  have  either  uniform  or  Poisson  distributions.  An  implementation  of  the  Mersenne 
twister  algorithm,  developed  by  Matsumoto  and  Nishimura  (i),  is  included.  The  output  from  the 
Mersenne  twister  is  used  to  generate  the  various  distributions  through  the  use  of  assorted 
transformation  algorithms. 

The  functions  presented  here  are  offered  as  compiler-independent  alternatives  to  functions 
defined  by  the  new  C++1 1  standard  (2).  Although  the  new  standard  defines  functions  for 
generating  pseudorandom  numbers  with  uniform  distributions,  normal  distributions,  Poisson 
distributions,  etc.,  it  does  not  specify  which  algorithms  will  be  used  to  generate  those 
pseudorandom  numbers.  Thus,  compilers  that  conform  to  the  new  standard  will  all  generate 
pseudorandom  numbers  from  the  various  distributions,  but  the  actual  numbers  that  are  generated 
may  differ  from  compiler  to  compiler. 

The  functions  presented  in  this  report  have  been  grouped  into  the  yRandom  namespace,  which  is 
summarized  at  the  end  of  this  report. 


2.  Generating  Pseudorandom  Integers  Using  the  Mersenne  Twister  19937 
Algorithm  -  The  Rand()  Function 


The  Rand()  function  uses  the  Mersenne  twister  19937  algorithm  to  generate  uniformly 
distributed  pseudorandom  integers  in  the  interval  [0,2  ).  According  to  Matsumoto  and 

1 QQ37 

Nishimura,  the  algorithm  has  a  period  of  2  -1  and  passes  the  Diehard  tests  for  statistical 

randomness. 

The  state  of  the  Mersenne  twister  is  stored  in  an  array  of  625  32-bit  unsigned  integers,  which  is 
passed  as  an  argument  to  the  Rand()  function.  The  initial  state  of  the  Mersenne  twister  can  be  set 
using  the  Initialize()  function  (see  section  3).  The  Initialize()  function  uses  a  user-supplied  32-bit 
integer  to  seed  a  simple  pseudorandom  number  generator  that  generates  the  initial  state  integers. 
This  effectively  allows  for  the  creation  of  multiple  pseudorandom  number  generators  (as  many 
as  2^^)  that  can  each  produce  independent,  reproducible  sequences  of  pseudorandom  integers. 
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The  code  contained  in  the  Rand()  function  differs  from  the  code  that  is  presented  by  Matsumoto 
and  Nishimura.  To  avoid  using  the  modulo  operator,  their  code  calculates  all  of  the  algorithm’s 
624  unsigned  state  integers  once  out  of  every  624  function  calls.  In  contrast,  the  Rand()  function 
uses  the  ternary  operator  to  avoid  using  the  modulo  operator. 

2.1  RandO  Code 


template<class  T>T  Rand(//<======================MERSENNE  TWISTER  (19937)  PRNG 

T  I[625]){//<-STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 
T  i=I [624] , j=i<623?i+l : 0,y=I [i]&0x80000000 | I [ j ]&0x7fffffff ; 
y=I[i]=I[i<227?i+397:i-227]''y>>l''(y&l)*0x9908b0clf,I[624]=j; 
return  y''(y''=(y''=(y''=y>>ll)<<7&0x9d2c5680)<<15&0xefc60000) >>18; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DAN2014 - 


2.2  RandO  Template  Classes 

T  T  should  be  a  32-bit  unsigned  integer  type.  Greater  than  32-bit  unsigned  integer 

types  can  be  used,  but  performance  may  suffer. 

2.3  RandO  Parameters 

I  I  points  to  a  625-element  array  that  contains  the  state  of  the  Mersenne  twister 

algorithm.  Each  time  the  RandO  function  is  called,  the  array  pointed  to  by  I  is 
modified.  The  initial  state  of  the  array  should  be  set  using  the  InitializeO  function 
(see  section  3). 

2.4  RandO  Return  Value 

'I'y 

The  RandO  function  returns  uniformly  distributed  pseudorandom  integers  in  the  interval  [0,2  ). 

2.5  RandO  Example 

The  following  example  uses  the  RandO  function  to  generate  one  billion  pseudorandom  integers. 

The  InitializeO  function,  introduced  in  section  3,  is  used  to  set  the  initial  state  of  the  Mersenne 

twister. 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_r'andom. h"// . yRandom 

int  main(){//<==========================EXAMPLE  FOR  THE  yRandom: : Rand ()  FUNCTION 


unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
for(int  i=l;i<1000000000;++i)yRandom: :Rand(I); 

printf("10*9  pseudorandom  numbers  generated  in  %.3f  seconds . \nThe  10*9th" 

"  number  is  %u. \n\n",double(clock())/CLOCKS_PER_SEC, yRandom: :Rand(I)); 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  AN  2  0 
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OUTPUT: 


10''9  pseudorandom  numbers  generated  in  1.778  seconds. 
The  number  is  2716480233. 


2.6  Comparison  of  the  Rand()  Function  to  Other  Mersenne  Twister  Implementations 


The  following  example  eompares  the  Rand()  funetion  to  Matsumoto  and  Nishimura’s  genrand() 
function,  as  well  as  to  the  C++1 1  built-in  Mersenne  twister  implementation.  The  example  code 
demonstrates  that  the  three  implementations  produce  identical  output  (at  least  for  the  first  10^ 
values  and  with  a  seed  value  of  1).  Note  that  time  values  will  vary  based  on  computer 
specifications,  compiler,  compiler  settings,  etc. 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  <random>// . MT19937 

#include  "matsumoto_nishimura . h"// . init_genrand( ) , genrand ( ) 

#include  "y_r'andom. h"// . yRandom 

int  main(){//<=====COMPARISON  BETWEEN  DIFFERENT  MERSENNE  TWISTER  IMPLEMENTATIONS 

double  t=0;// . elapsed  time 

init  genrand(l) : // . initialize  the  Matsumoto-Nishimura  implementation 

for (int  i=l; i< 1000000000; ++i) genrand ( ) ; 

printf( "Using  Matsumoto  and  Nishmimura's  genrand( ) : \n" ) ; 

printf("  10''9  pseudorandom  integers  generated  in  %.3f  seconds.  \n  The  10''9th" 
"  number  is  %u.\n\n", (clock()-t)/CLOCKS_PER_SEC,genrand()),t=clock(); 

std::mtl9937  g(l);// . initialize  the  C++11  built-in  implementation 

for (int  i=l; i<1000000000; ++i)g( ) ; 
printf ( "\nUsing  std: :mtl9937: \n"); 

printf("  10''9  pseudorandom  integers  generated  in  %.3f  seconds.  \n  The  10''9th" 
"  number  is  %u.\n\n", (clock()-t)/CLOCKS_PER_SEC,g()),t=clock(); 
unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//  init.  yRandom  implementation 
for(int  i=l;i<1000000000;++i)yRandom: :Rand(I); 
printf ( "\nUsing  yRandom: :Rand() : \n"); 

printf  ("  10''9  pseudorandom  integers  generated  in  %.3f  seconds.  \n  The  10''9th" 

"  number  is  %u.\n\n", (clock()-t)/CLOCKS_PER_SEC, yRandom: :Rand(I)); 
yRandom: : Initialize(I, 1) , init_genrand(l) ,g. seed(l) ; 
bool  check=true; 
for(int  i=0;i<1000000000;++i){ 
unsigned  x=yRandom: : Rand(I) ; 
if (x! =genrand( ) | | x! =g( ) )check=false; } 
printf ( "\nAre  the  first  10''9  pseudorandom  integers  generated  by\n 


"and  Nishimura's  genrand(),  std: :mtl9937, \n 
" ?  %s\n\n" , check? "YES" : "NO" ) ; 

}// - YAGENAUT@GMAI  L .  COM- 


Matsumoto 
and  yRandom: : Rand( )  identical'' 


LAST~UPDATED~28FEB2014' 
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OUTPUT: 


Using  Matsumoto  and  Nishmimura ' s  genrand(): 

10''9  pseudorandom  integers  generated  in  3.400  seconds. 
The  10''9th  number  is  2716480233. 


Using  std: :mtl9937: 

10''9  pseudorandom  integers  generated  in  2.231  seconds. 
The  10''9th  number  is  2716480233. 


Using  yRandom: :Rand() : 

10''9  pseudorandom  integers  generated  in  1.857  seconds. 
The  10''9th  number  is  2716480233. 


Are  the  first  10''9  pseudorandom  integers  generated  by 
Matsumoto  and  ISlishimura '  s  genrand(),  std:  :mtl9937, 
and  yRandom: : Rand( )  identical?  YES 


3.  Initializing  the  Mersenne  Twister  -  The  Initialize()  Function 


The  InitializeO  function  uses  an  algorithm  presented  by  Nishimura  and  Matsumoto  (3)  to 
initialize  the  625-element  array  that  is  used  to  store  the  state  of  the  Mersenne  twister  algorithm. 

3.1  InitializeO  Code 


templatecclass  T>void  Initialize(//<======INITIALIZE  STATE  OF  MERSENNE  TWISTER 

T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

unsigned  long  s){//< - SEED  [0,2''32) 

I [0] =s&0xffffffff , I [624] =0; 

for(int  i=l;i<624;++i)I[i]=(1812433253*(I[i-l]''I[i-l]>>30)+i)&0xffffffff; 

}//  /S//N//S//N/  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DAN2014 - 


3.2  InitializeO  Template  Classes 

T  T  should  be  a  32-bit  unsigned  integer  type.  Greater  than  32-bit  types  can  be  used, 

but  performance  may  suffer. 
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3.3  InitializeO  Parameters 


I  I  points  to  storage  for  a  625-element  array  that  is  used  to  store  the  state  of  the 

Mersenne  twister  algorithm. 

s  s  specifies  the  seed  for  the  algorithm  that  is  used  to  generate  the  initial  state  of  the 

TO 

Mersenne  twister  algorithm,  s  can  be  any  value  in  the  interval  [0,2  ). 


4.  Generating  Uniformly  Distributed  Pseudorandom  Numbers  -  The 
RandUO  Function 


The  RandUO  function  can  be  used  to  generate  uniformly  distributed  pseudorandom  numbers  that 
conform  to  the  probability  density  function  given  by  equation  1 . 


r  1 


/w  = 


b-a 


0, 


ioY a  <  x<b 
otherwise 


(1) 


Note  that  the  C-i-i-  standards  document  defines  a  slightly  different  probability  density  function  for 
uniformly  distributed  pseudorandom  numbers  (a<x<b  rather  than  a  <x<b).  The  decision  to 
exclude  the  lower  bound  from  the  distribution  rather  than  the  upper  bound  was  based  on  the 
belief  that  it  would  be  more  likely  to  help  the  user  to  avoid  a  divide-by-zero  error. 

The  RandUO  function  uses  the  transformation  presented  in  equation  2,  where  r  is  a  uniformly 
distributed  pseudorandom  number  in  the  interval  (0,1]. 

x  =  a  +  {b-  a)r  (2) 

Equation  3  can  be  used  to  generate  r  given  q,  a  uniformly  distributed  pseudorandom  integer  in 
the  interval  [0,2^^). 


4.1  RandUO  Code 


template<class  Txiouble  RandU(//<====UNIFORMLY  DISTRIBUTED  PSEUDORANDOM  DOUBLE 
T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

double  a, double  b){//< - LOWER  &  UPPER  BOUNDARIES  OF  DISTRIBUTION 

return  a+(b-a)*(Rand(I)+l.  )/4294967296;// . for  a=0  and  b=l,  (0,1] 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  AN  2  0 
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4.2  RandUO  Template  Classes 


T  T  should  be  a  32-bit  unsigned  integer  type.  Greater  than  32-bit  types  can  be  used, 

but  performance  may  suffer. 

4.3  RandUO  Parameters 

I  I  points  to  an  array  that  contains  the  state  of  the  Mersenne  twister  algorithm, 

a  a  specifies  a,  the  lower  bound  for  the  distribution, 

b  b  specifies  b,  the  upper  bound  for  the  distribution. 

4.4  RandUO  Return  Value 

The  RandUO  function  returns  uniformly  distributed  pseudorandom  numbers  in  the  interval  (a,b]. 
Note  that,  due  to  limitations  inherent  in  the  storage  of  double  precision  numbers,  the  return  value 
is  not  guaranteed  to  be  distinct  from  a  in  all  cases  (such  as  when  lal  »  la-bl). 

4.5  RandUO  Simple  Example 

The  following  example  uses  the  RandUO  function  to  generate  and  sum  one  billion  uniformly 
distributed  pseudorandom  numbers  in  the  interval  (2,6]. 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_random. h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom: :RandU()  FUNCTION 


unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandU(I,2,6); 
printf("10^9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f. \n\n",double(clock( ) )/CLOCKS_PER_SEC, s/1000000000) ; 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  AN  2  0 14*^ 


OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  8.255  seconds. 
Their  average  is  3.999977. 


4.6  RandUO  Binning  Example 

The  following  example  uses  the  RandU()  function  to  generate  and  bin  one  billion  uniformly 
distributed  pseudorandom  numbers  in  the  interval  (2,6]. 
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#include  <cstdio>// . printf() 

#include  "y_r'andom. h"// . yRandom 


int  main(){//<=================BINNING  EXAMPLE  FOR  THE  yRandom: :RandU()  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

int  M=1000000000; // . number  of  random  numbers  to  generate 

const  int  N=ll;// . number  of  bins  (not  counting  overflow  bin) 


double  B[N] ; /*<-*/for(int  i=0; i<N;++i)B[i]=4*double(i)/(N-l)+2; // . bins 

double  E[N+l]={0};/*<-*/for(int  i=l;i<N;++i)E[i]=M/(N-l. );//.. .expected  counts 

int  C[N+l];/*<-*/for(int  i=0;i<N+l;++i)C[i]=0;// . raw  counts 

for (int  i=0, j=0; i<M;++i,++C[ j ] , j=0) 

for(double  x=yRandom: :RandU(I,2,6);x>B[ j]&&j<N;++j); 
printfC  COUNT  COUNT\n  "); 

printf^'  BIN  ,  (RAW)  ,  (EXPECTED)  ,  %%DIFF  \n  "); 

printfC  - \n"); 

for(int  i=0;i<N+l;++i){ 

if(i==0)printf("  <=  %4.1f  ,",B[0]); 

else  if(i==N)printf("  >  %4.1f  ,",B[N-1]); 

else  printfC  (%4.1f  to  %4.1f]  , B[i-1] , B[i] ) ; 

printfC%lld  ,%13. If  " ,C[i] ,  E [i] ) ; 


E[i]>.l?printfC  ,%9 .4f%%\n",fabs(E [i] -C[i] )/E [i] *100) : printf C\n" ) ; } 
^11  E  N  AUT  ^)GI^AI  L  •  LAS  T^U  PD  AT  E  D^2  T  D  AN  2  0 


OUTPUT: 


COUNT 

COUNT 

BIN 

> 

(RAW) 

y 

(EXPECTED) 

y 

%DIFF 

<=  2 

0 

> 

0 

y 

0.0 

( 

2.0  to 

2.4] 

y 

99999019 

y 

100000000 . 0 

y 

0.0010% 

( 

2.4  to 

2.8] 

y 

100010759 

y 

100000000 . 0 

y 

0.0108% 

( 

2.8  to 

3.2] 

y 

99997646 

y 

100000000 . 0 

y 

0.0024% 

( 

3.2  to 

3.6] 

y 

100006913 

y 

100000000 . 0 

y 

0.0069% 

( 

3.6  to 

4.0] 

y 

100010417 

y 

100000000 . 0 

y 

0.0104% 

( 

4.0  to 

4.4] 

y 

99979386 

y 

100000000 . 0 

y 

0.0206% 

( 

4.4  to 

4.8] 

y 

99993880 

y 

100000000 . 0 

y 

0.0061% 

( 

4.8  to 

5.2] 

y 

99999652 

y 

100000000 . 0 

y 

0.0003% 

( 

5.2  to 

5.6] 

y 

100004658 

y 

100000000 . 0 

y 

0 . 0047% 

( 

5.6  to 

6.0] 

y 

99997670 

y 

100000000 . 0 

y 

0.0023% 

>  6 

0 

y 

0 

y 

0.0 
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5.  Generating  Normally  Distributed  Pseudorandom  Numbers  -  The  RandN() 
Function 


The  RandN()  function  can  be  used  to  generate  normally  distributed  pseudorandom  numbers  that 
conform  to  the  probability  density  function  given  by  equation  4,  where  jd  is  the  mean  of  the 
distribution  and  cr  is  the  standard  deviation. 

=  (4) 

ayl27r 

The  RandN()  function  uses  the  Box-Muller  transform  (4)  to  generate  normally  distributed 
pseudorandom  numbers: 

x  =  jU  +  <j.J-  21n(r^)  cos(2;zr^ )  (5) 

where  and  are  uniformly  distributed  pseudorandom  numbers  in  the  interval  (0,1].  Equation  3 
can  be  used  to  generate  and  n  given  qa  and  qt,  two  uniformly  distributed  pseudorandom 
integers  in  the  interval  [0,2  ). 

5.1  RandNO  Code 

template<class  T>double  RandN(//<=====ISIORMALLY  DISTRIBUTED  PSEUDORANDOM  DOUBLE 
T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

double  m, double  s){//< - MEAN  &  STANDARD  DEVIATION  OF  DISTRIBUTION 

return  m+s*sqrt(-2*log((Rand(I)+l. )/4294967296) ) 
*cos(1.4629180792671596E-9*(Rand(I)+l.)); 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1  D  AN  2  9 


5.2  RandNO  Template  Classes 

T  T  should  be  a  32-bit  unsigned  integer  type.  Greater  than  32-bit  types  can  be  used, 

but  performance  may  suffer. 

5.3  RandNO  Parameters 

I  I  points  to  an  array  that  contains  the  state  of  the  Mersenne  twister  algorithm, 

m  m  specifies  /l,  the  mean  of  the  distribution, 

s  s  specifies  cr,  the  standard  deviation  of  the  distribution. 
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5.4  RandNO  Return  Value 

The  RandNO  function  returns  normally  distributed  pseudorandom  numbers. 

5.5  RandNO  Simple  Example 

The  following  example  uses  the  RandNO  function  to  generate  and  sum  one  billion  normally 
distributed  pseudorandom  numbers  with  fu  =  4  and  a  =  0.5. 

#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_r'andom. h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom: ;RandN()  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandN(I,4, .5); 
printf ( "10''9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f. \n\n",double(clock( ) )/CLOCKS_PER_SEC, s/1000000000) ; 

^11  E  N  AUT  ^)GI^AI  L  •  LAS  T^U  PD  AT  E  D^2  T  D  AN  20 14*'^*"*’*"*’*"^*^“^ 


OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  46.213  seconds. 
Their  average  is  3.999998. 


5.6  RandNO  Binning  Example 

The  following  example  uses  the  RandNO  function  to  generate  and  bin  one  billion  normally 
distributed  pseudorandom  numbers  with  ju  =  4  and  cr  =  0.5.  The  ErfO  function  is  an 
implementation  of  an  algorithm  presented  by  Abramowitz  and  Stegun  (5). 
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#include  <cstdio>// . printf() 

#include  "y_r'andom. h"// . yRandom 


inline  double  Erf (//<============================RETURI\IS  THE  ERROR  FUNCTION  OF  X 

double  x){//< - ANY  REAL  NUMBER 

double  t=l/(l+.3275911*fabs(x));//. .see  Abramowitz  and  Stegun,  7.1.26  (p.  299) 
double  a [ ] ={ . 254829592, - . 284496736, 1.421413741, -1.453152027,1 . 061405429}; 
double  s=0;/*<-*/for(int  i=0;i<5;++i)s+=a[i]*pow(t,i+l); 

return  (x<0?-l:l)*(l-s*exp(-x*x));// . note  erf ( -x)  =  -erf (x) 

^11  E  N  AUT  ^)GI^AI  L  •  LAS  T  PDAT  E  D^2 1 D  AN  2  0 

int  main(){//<=================BINNING  EXAMPLE  FOR  THE  yRandom; ;RandN()  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

int  M=1000000000; // . number  of  random  numbers  to  generate 

const  int  N=ll;// . number  of  bins  (not  counting  overflow  bin) 

double  B[N] ; /*<-*/for(int  i=0; i<N;++i)B[i]=4*double(i)/(N-l)+2; // . bins 

double  E[N+l]={0};/*<-*/E[0]=E[N]=.5*(l+Erf((B[0]-4)/sqrt(2*.5*.5)))*M;//.exp. 
for (int  i=l; i<N;++i)E [i]=( . 5*(1+Erf ( (B[i] -4)/ 

sqrt(2*.5*.5)))-.5*(l+Erf((B[i-l]-4)/sqrt(2*.5*.5))))*M; 

int  C[N+l];/*<-*/for(int  i=0;i<N+l;++i)C[i]=0;// . raw  counts 

for (int  i=0, j=0; i<M;++i,++C[ j ] , j=0) 

for(double  x=yRandom: :RandN(I,4, .5);x>B[ j]&&j<N;++j); 


printfC  COUNT  COUNT\n  "); 

printf^'  BIN  ,  (RAW)  ,  (EXPECTED)  ,  %%DIFF  \n  "); 

printf("  - \n"); 

for(int  i=0;i<N+l;++i){ 

if(i==0)printf("  <=  %4.1f  ,",B[0]); 

else  if (i==N)printf ( "  >  %4.1f  ,",B[N-1]); 

else  printfC  (%4.1f  to  %4.1f]  , " , B[i-1] , B[i] ) ; 


printf("%lld  ,%13. If " ,C[i] , E [i] ) ; 

E[i]>.l?printf("  ,%9 .4f%%\n",fabs(E [i] -C[i] )/E [i] *100) : printf ( "\n" ) ; } 

^11  YAG  E  N  AUT  ^^GI^AI  L  •  LAS  T^U  PDAT  E  D'^  2  T  D  AN  2  0 1 4^*^^*^*^*^ 


OUTPUT: 


COUNT 

COUNT 

BIN 

> 

(RAW) 

y 

(EXPECTED) 

y 

%DIFF 

<=  2 

0 

y 

31839 

y 

31686.0 

y 

0.4827% 

( 

2.0  to 

2.4] 

y 

655830 

y 

655516.1 

y 

0.0479% 

( 

2.4  to 

2.8] 

y 

7507620 

y 

7510327.1 

y 

0.0360% 

( 

2.8  to 

3.2] 

y 

46597020 

y 

46601762.1 

y 

0.0102% 

( 

3.2  to 

3.6] 

y 

157066112 

y 

157056047.3 

y 

0 . 0064% 

( 

3.6  to 

4.0] 

y 

288150018 

y 

288144661.9 

y 

0.0019% 

( 

4.0  to 

4.4] 

y 

288145262 

y 

288144660.9 

y 

0.0002% 

( 

4.4  to 

4.8] 

y 

157043835 

y 

157056047.3 

y 

0.0078% 

( 

4.8  to 

5.2] 

y 

46607503 

y 

46601762.1 

y 

0.0123% 

( 

5.2  to 

5.6] 

y 

7508108 

y 

7510327.1 

y 

0.0295% 

( 

5.6  to 

6.0] 

y 

655295 

y 

655516.1 

y 

0.0337% 

>  6 

0 

y 

31558 

y 

31686.0 

y 

0.4041% 
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6.  Generating  Uniformly  Distributed  Pseudorandom  Integers  -  The  Randl() 
Function 


The  Randl()  function  can  be  used  to  generate  uniformly  distributed  pseudorandom  integers  that 
conform  to  the  probability  density  function  given  by  equation  6. 


f{k)  = 


- ^ - ,  for  a<k<b 

b-a  +  l 

0,  otherwise 


(6) 


The  Randl()  function  uses  the  transformation  presented  in  equation  7,  where  q  is  &  uniformly 


distributed  pseudorandom  integer  in  the  interval  [0,2  ). 


k  =  a  +  trunc 


^{b-a  +  l) 


(7) 


6.1  RandIO  Code 


template<class  T>long  RandI(//<=========UNIFORMLY  DISTRIBUTED  PSEUDORANDOM  INT 

T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

long  a, long  b){//< - LOWER  &  UPPER  BOUNDARIES  OF  DISTRIBUTION 

return  a+T(Rancl(I)/4294967296.*(b-a+l));// . [a^b] 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DAN2014 - 


6.2  RandIO  Template  Classes 

T  T  should  be  a  32-bit  unsigned  integer  type.  Greater  than  32-bit  types  can  be  used, 

but  performance  may  suffer. 

6.3  RandIO  Parameters 

I  I  points  to  an  array  that  contains  the  state  of  the  Mersenne  twister  algorithm, 

a  a  specifies  a,  the  lower  bound  for  the  distribution, 

b  b  specifies  b,  the  upper  bound  for  the  distribution. 

6.4  RandIO  Return  Value 

The  RandIO  function  returns  a  uniformly  distributed  pseudorandom  integer  in  the  interval  [a,b] . 

6.5  RandIO  Simple  Example 

The  following  example  uses  the  Randl()  function  to  generate  and  sum  one  billion  uniformly 
distributed  pseudorandom  integers  in  the  interval  [-5,4]. 
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#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_random. h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom: :RandI()  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandI(I, -5,4); 
printf ( "10''9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f. \n\n",double(clock( ) )/CLOCKS_PER_SEC, s/1000000000) ; 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  AN  2  0 


OUTPUT: 


10^9  pseudorandom  numbers  generated  and  summed  in  10.140  seconds. 
Their  average  is  -0.500057. 


6.6  RandlO  Binning  Example 

The  following  example  uses  the  Randl()  function  to  generate  and  bin  one  billion  uniformly 
distributed  pseudorandom  integers  in  the  interval  [-5,4]. 


#include  <cstdio>// . printf() 

#include  "y_random. h"// . yRandom 


int  main(){//<=================BINNING  EXAMPLE  FOR  THE  yRandom: :RandI()  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

int  M=1000000000; // . number  of  random  numbers  to  generate 

const  int  N=ll;// . number  of  bins  (not  counting  overflow  bin) 


double  B[N];/*<-*/for(int  i=0;i<N;++i)B[i]=i-6;// . bins 

double  E[N+l]={0};/*<-*/for(int  i=l;i<N;++i)E[i]=M/(N-l. );//.. .expected  counts 

int  C[N+l];/*<-*/for(int  i=0;i<N+l;++i)C[i]=0;// . raw  counts 

for (int  i=0, j=0; i<M;++i,++C[ j ] , j=0) 

for(double  x=yRandom: :RandI(I, -5,4) ;x>B[ j ]&&j<N;++j ) ; 
printfC  COUNT  COUNT\n  "); 

printf^'  BIN  ,  (RAW)  ,  (EXPECTED)  ,  %%DIFF  \n  "); 

printfC  - \n"); 

for(int  i=0;i<N+l;++i){ 

if (i==0)printf ( "  <=  %4.1f  ,",B[0]); 

else  if(i==N)printf("  >  %4.1f  ,",B[N-1]); 

else  printfC  (%4.1f  to  %4.1f]  , " , B[i-1] , B[i] ) ; 

printfC%lld  ,%13. If  " ,C[i] ,  E [i] ) ; 


E[i]>.l?printfC  ,%9 .4f%%\n",fabs(E [i] -C[i] )/E [i] *100) : printf C\n" ) ; } 

^11  E  N  AUT  ^)GI^AI  L  •  LAS  T  ^U  PDAT  E  D^2  T  D  AN  2  0 14^^*^^ 
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OUTPUT: 


COUNT 

COUNT 

BIN 

> 

(RAW) 

> 

(EXPECTED) 

> 

%DIFF 

<=  -6.0 

> 

0 

> 

0.0 

(-6.0  to  -5.0] 

> 

99999019 

> 

100000000 . 0 

> 

0.0010% 

(-5.0  to  -4.0] 

> 

100010759 

> 

100000000 . 0 

> 

0.0108% 

(-4.0  to  -3.0] 

> 

99997646 

> 

100000000 . 0 

> 

0.0024% 

(-3.0  to  -2.0] 

> 

100006913 

> 

100000000 . 0 

> 

0.0069% 

(-2.0  to  -1.0] 

> 

100010417 

> 

100000000 . 0 

> 

0.0104% 

(-1.0  to  0.0] 

> 

99979386 

> 

100000000 . 0 

> 

0.0206% 

(  0.0  to  1.0] 

> 

99993881 

> 

100000000 . 0 

> 

0.0061% 

(  1.0  to  2.0] 

> 

99999651 

> 

100000000 . 0 

> 

0.0003% 

(  2.0  to  3.0] 

> 

100004658 

> 

100000000 . 0 

> 

0 . 0047% 

(  3.0  to  4.0] 

> 

99997670 

> 

100000000 . 0 

> 

0.0023% 

>  4.0 

> 

0 

> 

0.0 

7.  Generating  Poisson-Distributed  Pseudorandom  Integers  -  The  RandP() 
Function 


The  RandP()  function  can  be  used  to  generate  Poisson  distributed  pseudorandom  integers  that 
conform  to  the  probability  density  function  given  by  equation  8,  where  jj  is  the  mean  of  the 
distribution. 

fik)  =  (8) 

k\ 

The  RandP()  function  uses  a  transformation  described  by  Knuth  (6)  to  generate  Poisson- 
distributed  pseudorandom  integers  by  finding  the  largest  k  that  satisfies  equation  9. 

(9) 

1=0 

where  r,  is  a  uniformly  distributed  pseudorandom  number  in  the  interval  (0,1].  Equation  3  can  be 
used  to  generate  each  r,  given  qi,  a  uniformly  distributed  pseudorandom  integer  in  the  interval 
[0,2'"). 

Note  that  the  for  large  /r,  the  RandP()  function  is  slow.  This  problem  can  be  overcome  by 
making  use  of  the  fact  that,  for  large  /J,  Poisson  distributions  are  approximately  normal. 
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7.1  RandPO  Code 


template<class  T>T  RandP(//<=====POISSON-DISTRIBUTED  PSEUDORANDOM  UNSIGNED  INT 
T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

double  m){//< - MEAN  (MUST  BE  GREATER  THAN  ZERO) 

T  k=0;/*<-*/for(double  P=l,E=exp(-m);P>E;P*=(Rand(I)+l. )/4294967296)++k; 
return  k-1; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21IAN2014 - 


7.2  RandPO  Template  Classes 

T  T  should  be  a  32-bit  unsigned  integer  type.  Greater  than  32-bit  types  can  be  used, 

but  performance  may  suffer. 

7.3  RandPO  Parameters 

I  I  points  to  an  array  that  contains  the  state  of  the  Mersenne  twister  algorithm, 

m  m  specifies  jd,  the  mean  of  a  Poisson  distribution. 

7.4  RandPO  Return  Value 

The  RandPO  function  returns  a  Poisson-distributed  pseudorandom  unsigned  integer. 

7.5  RandPO  Simple  Example 

The  following  example  uses  the  RandPO  function  to  generate  and  sum  one  billion  Poisson- 
distributed  pseudorandom  integers  with  =  1 . 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_r'andom. h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom: :RandP()  FUNCTION 


unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandP(I,l); 
printf("10*9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f. \n\n",double(clock( ) )/CLOCKS_PER_SEC, s/1000000000) ; 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  D^2 1 D  AN  2  0 


OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  35.708  seconds. 
Their  average  is  1.000013. 
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7.6  RandPO  Binning  Example 


The  following  example  uses  the  RandP()  function  to  generate  and  bin  one  billion  Poisson 
distributed  pseudorandom  integers  with  [A  =  \. 


#include  <cstdio>// . printf() 

#include  "y_random. h"// . yRandom 


int  main(){//<=================BINNING  EXAMPLE  FOR  THE  yRandom: :RandP()  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

int  M=1000000000; // . number  of  random  numbers  to  generate 

const  int  N=ll;// . number  of  bins  (not  counting  overflow  bin) 

double  B[N];/*<-*/for(int  i=0;i<N;++i)B[i]=i;// . bins 

double  E[N+1]={0};// . expected  counts 

double  F=l;/*<-*/for(int  k=0; l«N+l;++k, F*=k)E[k]=pow(l. , k)*exp( -1. )/F*M; 

int  C[N+l];/*<-*/for(int  i=0;i<N+l;++i)C[i]=0;// . raw  counts 

for (int  i=0, j=0; i<M;++i,++C[ j ] , j=0) 

for(double  x=yRandom: :RandP(I,l);x>B[ j]&&j<N;++j); 


printfC  COUNT  COUNT\n  "); 

printf^'  BIN  ,  (RAW)  ,  (EXPECTED)  ,  %%DIFF  \n  "); 

printfC'  - \n"); 

for(int  i=0;i<N+l;++i){ 

if(i==0)printf("  <=  %4.1f  C',B[0]); 

else  if(i==N)printf("  >  %4.1f  ,",B[N-1]); 

else  printfC  (%4.1f  to  %4.1f]  , B[i-1] , B[i] ) ; 


printf("%lld  ,%13. If " ,C[i] > E [i] ) ; 

E[i]>.l?printf("  ,%9 .4f%%\n",fabs(E [i] -C[i] )/E [i] *100) : printf ( "\n" ) ; } 
^11  YAG  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  2 1 D  AN  2  0 


OUTPUT: 


COUNT 

COUNT 

BIN 

y 

(RAW) 

y 

(EXPECTED) 

y 

%DIFF 

<=  0 

.0 

y 

367886090 

y 

367879441.2 

y 

0.0018% 

( 

0.0  to 

1.0] 

y 

367871283 

y 

367879441.2 

y 

0.0022% 

( 

1.0  to 

2.0] 

y 

183927644 

y 

183939720.6 

y 

0.0066% 

( 

2.0  to 

3.0] 

y 

61321837 

y 

61313240.2 

y 

0.0140% 

( 

3.0  to 

4.0] 

y 

15332140 

y 

15328310.0 

y 

0.0250% 

( 

4.0  to 

5.0] 

y 

3067729 

y 

3065662.0 

y 

0.0674% 

( 

5.0  to 

6.0] 

y 

510230 

y 

510943.7 

y 

0.1397% 

( 

6.0  to 

7.0] 

y 

72871 

y 

72992.0 

y 

0.1657% 

( 

7.0  to 

8.0] 

y 

9053 

y 

9124.0 

y 

0.7781% 

( 

8.0  to 

9.0] 

y 

999 

y 

1013.8 

y 

1.4576% 

( 

9.0  to 

10.0] 

y 

112 

y 

101.4 

y 

10.4779% 

>  10 

.0 

y 

12 

y 

9.2 

y 

30.2061% 
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8.  Code  Summary 


A  summary  sheet  is  provided  at  the  end  of  this  report.  It  presents  the  yRandom  namespace, 
which  contains  the  six  functions  that  are  described  in  this  report. 
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yRandom  Summary 


y_random. h 


Rand()  Example 


#ifndef  Y_RANDOM_GUARD//  See  YagePj  R.3.  "Generating  Psuedorandom  Numbers 

#define  Y_RANDOM_GUARD//  from  Various  Distributions  Using  C++" 

#include  <cmath>// . exp( ) j log( ) j sin( ) j  sqrt() 

template<class  T>void  InitRFH?S(//<======INITIALIZE  STATE  OF  MERSENNE  TWISTER 

T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

unsigned  long  s){//< - SEED  [0,2''32) 

I [ 0 ] =s&0xf f f f f f f f J I [ 624 ] =0; 

for(int  i=l;i<624;++i)I[i]=(1812433253*(I[i-l]''I[i-l]>>30)+i)&0xffffffff; 

}//~~~~YAGENAUT@GMAI1^C0M'- - - - - - - - - ~~LAST~UPDATED~21:AN2014— - 

template<class  T>T  Rand(//<======================MERSENNE  TWISTER  (19937)  PRNG 

T  I[625]){//<-STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 
T  i=I[624], j=i<623?i+l:0,y=I[i]&0x80000000|l[j]&0x7fffffff; 
y=I[i]=I[i<227?i+397:i-227]''y>>l''(y&l)*0x9908b0dfjI[624]=j; 
return  y'(y''=(y''=(y''=y>>ll)<<7&0x9d2c5680)<<15&0xefc60000)>>18; 

}//~^^~YAGENAUT^MAIL.C0M - - - ^^-.^^^^^~^^^^^lAST~UPDATED~213AN2014— - - 

template<class  T>double  RandU(//<====UNIFORMLY  DISTRIBUTED  PSEUDORANDOM  DOUBLE 
T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

double  a, double  b){//< . LOWER  &  UPPER  BOUNDARIES  OF  DISTRIBUTION 

return  a+(b-a)*(Rand(I)+l.  )/4294967296j // . for  a=0  and  b=l,  (0jl] 

}//~^^~YAGENAUT@GMAIL.C0M~'— - - - - - LAST~UPDATED~213AN2014 - 

template<class  T>double  RandN(//<«===NORMALLY  DISTRIBUTED  PSEUDORANDOM  DOUBLE 
T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

double  m, double  s){//< . MEAN  &  STANDARD  DEVIATION  OF  DISTRIBUTION 

return  m+s*sqrt(-2*log( (Rand(I)+l . ) 74294967296)) 

*cos(1.4629180792671596E-9*(Rand(I)+l.)); 

}//~~~~YAGENAUT@GMAIL.C0M - LAST~UPDATED~21:AN2014~~~~'-'- 

template<class  T>long  RandI(//<=========UNIFORMLY  DISTRIBUTED  PSEUDORANDOM  INT 

T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

long  a, long  b){//< . LOWER  &  UPPER  BOUNDARIES  OF  DISTRIBUTION 

return  a+T(Rand(I)/4294967296.*(b-a+l));// . [a^b] 

}// - YAGENAUT@GMAIL.COM - ~'-~LAST-UPDATED~21:AN2014 - 

template<class  T>T  RandP(//<==«=POISSON-OISTRIBUTED  PSEUDORANDOM  UNSIGNED  INT 
T  I[625],//<--STATE  OF  MERSENNE  TWISTER  (FOR  T,  USE  A  32-BIT  UNSIGNED  INT) 

double  m){//< . MEAN  (MUST  BE  GREATER  THAN  ZERO) 

T  k=0;/*<-*/for (double  P=l, E=exp( -m) jP>E;P*=(Rand(I)+l. )/4294967296)++k; 
return  k-1; 

}//— YAGENAUT@GMAIL  .  COM - - - LAST'-UPDATED'-21:AN2014 - 


#endif 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) jCLOCKS_PER_SEC 

#include  "y_random.  h"// . yRandom 

int  main(){//<====««=============«===EXAMPLE  FOR  THE  yRandom:  :Rand()  FUNCTION 


unsigned  I[625]j/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
for (int  i=lj i< 1000000000; ++i)yRandom: :Rand(I); 

printf  (”10''9  pseudorandom  numbers  generated  in  %.3f  seconds.  \nThe  10''9th" 

"  number  is  %u.\n\n"jdouble(clock())/CLOCKS_PER_SECjyRandom: :Rand(I)); 

}//..— - YAGENAUT@GMAIL.COM~- - - - LAST~UPDATED~21:AN2014 - 

OUTPUT: 


10''9  pseudorandom  numbers  generated  in  1.778  seconds. 
The  10''9th  number  is  2716480233. 


RandU()  Probability  Density  Function  and  Simple  Example 


,  — ,  for  a  <  x<b 

f{x)  =  \h-a 

0,  otherwise 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) j  CLOCKS_PER_SEC 

#include  ”y_random.  h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom: : RandU( )  FUNCTION 


unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandU(I,2j6); 
printf ("10''9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f.\n\n",double(clock())/CLOCKS_PER_SEC, s/1000000000); 

}// - ^yAGENAUT@GMAIL.COM~~~~~~~~~~~~~~ - - - LAST~UPDATED~21]AN2014~~ - 


OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  8.255  seconds. 
Their  average  is  3.999977. 


RandN()  Probability  Density  Function  and  Simple  Example 


RandN()  Binning  Example 


#include  <cstdio>// . printf() 

#in elude  "y_random. h"// . yRandom 


inline  double  Erf (//<============================RETURNS  THE  ERROR  FUNCTION  OF  X 

double  x){//< . ANY  REAL  NUMBER 

double  t=l/(l+.3275911*fabs(x));//. .see  Abramowitz  and  Stegunj  7.1.26  (p.  299) 
double  a[]={.254829592j -.284496736^1.421413741^-1.453152027,1.061405429}; 
double  s=0;/*<-*/for(int  i=0;i<5;++i)s+=a[i]*pow(t,i+l); 

return  (x<0?-l:l)*(l-s*exp(-x*x));// . note  erf (-x)  =  -erf (x) 

}// - —YAGENAUT@GMAIL.COM - - - LAST~UPDATED~21OAN2014~~~~~~ 

int  main(){//<=================BINNING  EXAMPLE  FOR  THE  yRandom: : RandN( )  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

int  M=1000000000; // . number  of  random  numbers  to  generate 

const  int  N=ll;// . number  of  bins  (not  counting  overflow  bin) 


double  B[N];/*<-*/for(int  i=0;i<N;++i)B[i]=4*double(i)/(N-l)+2;// . bins 

double  E[N+l]={0};/*<-*/E[0]=E[N]=.5*(l+Erf((B[0]-4)/sqrt(2».5*.5)))*M;//.exp. 
for (int  i=l;i<N;++i)E[i]=(.5*(l+Erf((B[i]-4)/ 

sqrt(2*.5».5)))-.5*(l+Erf((B[i-l]-4)/sqrt(2*.5*.5))))*M; 

int  C[N+l];/*<-*/for(int  i=0;i<N+l; ++i)C[i]=0; // . raw  counts 

for (int  i=0j j=0;i<M;++ij++C[ j], j=0) 

for(double  x=yRandom: :RandN(I,4, .5);x>B[j]&&j<N;++j); 
printf ("  COUNT  COUNT\n  "); 

printff'  BIN  ,  (RAW)  ,  (EXPECTED)  ,  %%DIFF  \n  "); 

printf("  - \n"); 

for(int  i=0;i<N+l;++i){ 

if(i==0)printf ("  <=  %4.1f  j"jB[0]); 

else  if(i==N)printf("  >  %4.1f  ,",B[N-1]); 

else  printf("  (%4.1f  to  %4.1f]  , B[i-1] ,B[i] ) ; 

printf ("%lld  ,%13.1f",C[i],E[i]); 


E[i]>.l?printf("  ,%9.4f%%\n",fabs(E[i] -C[i] )/E[i] *100) : printf (" \n"); } 

}// - —YAGENAUT@GMAIL.COM - - - - - - - LAST~UPDATED~21:AN2014 - 

OUTPUT: 


/W  = 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ),CLOCKS_PER_SEC 

#include  "y_random.  h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom :: RandN( )  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandN(I,4j .5); 
printf (”10''9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f.\n\n",double(clock())/CLOCKS_PER_SEC, s/1000000000); 
}//.^— ~^YAGENAUT@GMAIL.C0M~~ - - - - - LAST~UPDATED~21:AN2014 - 


OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  46.213  seconds. 
Their  average  is  3.999998. 


Randl()  Probability  Density  Function  and  Simple  Example 


1 


f(k)  =  >.b-a^\ 

[  0. 


for  «  <  ^  ^ 
otherwise 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ),CLOCKS_PER_SEC 

#include  "y_random.  h"// . yRandom 

int  main(){//<==================SIMPLE  EXAMPLE  FOR  THE  yRandom :: Randl( )  FUNCTION 


unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandI(I, -5,4); 
printf ("10''9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f.\n\n",double(clock())/CLOCKS_PER_SEC, s/1000000000); 

}// - ~YAGENAUT@GMAIL.COM~~~~~~~~~~~~~~ - - - ~LAST~UPDATED~21:AN2014 - ~~~ 


BIN 

COUNT 

,  (RAW)  , 

COUNT 

(EXPECTED)  , 

%DIFF 

<=  2 

•0 

,  31839  , 

31686.0  , 

0.4827% 

2.0  to 

2.4]  , 

,  655830  , 

655516.1  , 

0.0479% 

2.4  to 

2.8]  , 

,  7507620  , 

7510327.1  , 

0.0360% 

2.8  to 

3.2]  , 

,  46597020  , 

46601762.1  , 

0.0102% 

3.2  to 

3.6]  , 

,  157066112  , 

,  157056047.3  , 

0.0064% 

3.6  to 

4.0]  , 

,  288150018  , 

,  288144661.9  , 

0.0019% 

4.0  to 

4.4]  , 

,  288145262  , 

,  288144660.9  , 

0.0002% 

4.4  to 

4.8]  , 

,  157043835  , 

,  157056047.3  , 

0.0078% 

4.8  to 

5.2]  , 

,  46607503  , 

46601762.1  , 

0.0123% 

5.2  to 

5.6]  , 

,  7508108  , 

7510327.1  , 

0.0295% 

5.6  to 

6.0]  J 

,  655295  , 

655516.1  , 

0.0337% 

>  6 

.0 

,  31558  , 

31686.0  , 

0.4041% 

OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  10.140  seconds. 
Their  average  is  -0.500057. 


RandP()  Probability  Density  Function  and  Simple  Example 


f{k)  = 


H  e 


k\ 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ),CLOCKS_PER_SEC 

#include  ”y_random.  h"// . yRandom 

int  main(){//<=*==«=*=====«===SIMPLE  EXAMPLE  FOR  THE  yRandom :: RandP( )  FUNCTION 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double  s=0;/*<-*/for(int  i=0;i<1000000000;++i)s+=yRandom: :RandP(I,l); 
printf ("10''9  pseudorandom  numbers  generated  and  summed  in  %.3f  seconds. \n" 
"Their  average  is  %f.\n\n",double(clock())/CLOCKS_PER_SEC, s/1000000000); 

}// - YAGENAUT@GMAIL.COM~~ - - - - - LAST~UPDATED~21:AN2014 - 


OUTPUT: 


10''9  pseudorandom  numbers  generated  and  summed  in  35.708  seconds. 
Their  average  is  1.000013. 


17 


9.  References 


1.  Matsumoto,  M.;  Nishimura,  T.  Mersenne  Twister:  A  623 -Dimensionally  Equidistributed 
Uniform  Pseudo-Random  Number  Generator.  ACM  Transactions  on  Modeling  and 
Computer  Simulation  1998,  8,  (1),  3-20. 

2.  American  National  Standards  Institute  Standard  INCTTS/ISO/IEC  14882-2011.  Information 
Technology  -  Programming  Languages  -  C++  2012. 

3.  Nishimura,  T;  Matsumoto,  M.  A  C-Program  for  MT19937,  With  Initialization  Improved, 

26  January  2002.  http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES 
/mtI9937ar.c.  (accessed  12  December  2013). 

4.  Box,  G.  E.;  Muller,  M.  E.  A  Note  on  the  Generation  of  Random  Normal  Deviates.  The 
Annals  of  Mathematical  Statistics  1958,  29  (2),  610-611. 

5.  Abramowitz,  M.,  Stegun,  I.  A.,  Eds.  Handbook  of  Mathematical  Functions  With  Formulas, 
Graphs,  and  Mathematical  Tables;  Applied  Mathematics  Series  55;  U.S.  National  Bureau  of 
Standards:  Washington,  DC,  December  1972,  section  7.1.26. 

6.  Knuth,  D.  E.  Seminumerical  Algorithms,  the  Art  of  Computer  Programming;  Vol.  2; 

Addison  Wesley:  Reading,  MA,  1997. 


18 


1  DEFENSE  TECHNICAL 
(PDF)  INFORMATION  CTR 

DTIC  OCA 

2  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
RDRL  CIO  LL 

IMAL  HRA  MAIL  &  RECORDS  MGMT 

1  GOVT  PRINTG  OFC 
(PDF)  AMALHOTRA 

1  RDRL  WML  A 
(PDF)  R  YAGER 


19 


Intentionally  left  blank. 


20 


